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EFFECTS OF CAVITATION ON UNDERWATER 
SHOCK LOADINGS - Part I 

Abstract 

Reported here are analytic formulations, together with 
one-dimensional results, in an investigation of the title 
subject. It is shown that either displacement or a displacement 
potential may be used as the basic dependent variable for a 
finite element analysis. Artificial damping is found to be 
needed to suppress spurious oscillations (a numerical phenom- 
enon) near cavity boundaries. Adequacy of the method is demon- 
strated by comparison with published results of Bleich and 
Sandler. Some results are given for effects of cavitation on 
the performance of resilient attenuators. 
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EFFECTS OF CAVITATION ON UNDERWATER 
SHOCK LOADINGS - Part I 

1. Introduction 

1.1 Earlier Work . Motivation for the present investigation origin- 
ated with the "shock shield" proposed by Geers in Ref. 1. The 
shield is a gas-filled cushion (GFC) to be fitted to the exterior 
of a submarine hull. If the cushion has sufficient thickness 

it can greatly reduce the magnitude of underwater shock loads 
transmitted to the hull. Related concepts involving the applica- 
tion of resilient elastic layers (REL) are treated by Geers in 
Ref. 2. 

1.2 Cavitation Effects . The two-dimensional analyses of Refs. 1 
and 2 neglect possible effects of cavitation. It is well-known 
(see Refs. 3 and 4), however, that a highly compliant submerged 
object will produce a negative pressure scattered wave in response 
to an incident shock wave. If the shock pressure is much greater 
than the hydrostatic pressure, cavitation will be induced in the 
fluid. Such cavitation may significantly increase the shock load- 
ing on the submerged body. It is thus evident that an adequate 
investigation of the effectiveness of resilient attenuators 
requires evaluation of effects of cavitation on performance. 

1.3 Plan for Investigation . The present investigation is 
divided into two phases. The first phase is the subject of this 
report. It involves consideration of representative one- 
dimensional problems for the purpose of determining the relative 
merits of alternate choices for: dependent variables, time 
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integration algorithms, and spatial and temporal discretization. 
The one -dimensional context allows rapid and inexpensive compu- 
tations and allows comparison of results with those reported by 
others. An account of this phase is given in the remaining 
sections of this report. 

The second phase of the investigation consists of extensions 
to problems in three dimensions. Reasonable limitations on com- 
puter core capacity and processing time require that the problems 
be axisymmetric and, thus, mathematically two-dimensional. 

1.4 Fluid Model . It is known that fluids do have some capacity 
for sustaining negative pressure (tension) . Some data are given 
in Ref. 5. The influence of dissolved gas on the development 
of cavitation is considered in Ref. 6. For the purpose of the 
present investigation it is advantageous, and presumably con- 
servative, to assume that the transition from the normal to the 
cavitated state takes place without delay when the absolute 
pressure reaches zero. 

In the initial stages of this investigation the fluid was 
treated as bilinear with a greatly reduced bulk modulus in the 
negative pressure region. Subsequent developments disclosed 
that the expected advantages of the bilinear model were not 
achieved and the bulk modulus was henceforth assumed to be zero 
in the cavitated region. 



4 



2. Choice of Dependent Variable 

2.1 Failure of the Pressure Formulation . At the outset, this 
investigator expected that a formulation of the governing 
equations using fluid pressure p as the basic dependent variable 
would be advantageous. This expectation was based on previous 
successful finite element applications to propagation problems 
(e.g., see Refs. 7-9). Prior applications did not involve cavi- 
tation, but the bilinear fluid model was expected to handle 
successfully cavitation effects. 

At an early stage of the investigation, duplication of the 
results of the example problem of Ref. 4 was attempted. These 
trials gave solutions which correctly tracked the growth of the 
cavitated region, but failed to show its subsequent contraction 
and collapse. Efforts to discover the reason for the failure of 
the pressure formulation led to a simple test problem which 
determines whether a proposed formulation can correctly track 
the contraction of a cavitated region. The problem is defined 
in the following section. 

2.2 "Water-Hammer" Problem. The rapid pressure rise which 
accompanies the sudden interruption of water flow in a closed 
conduit is known as water-hammer. We here consider a flow in a 
zero pressure (or a small negative pressure, if the bilinear 
fluid model is used) cavitated region with positive dilatation 

e . Thus we have, for a semi - inf ini te region x > 0, the initial 
values : 

p(x,0) = 0, (pressure) 
e(x,0) = e Q > 0, (dilatation) 

u(x,0) = -v q < 0. (velocity) 
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The boundary condition at x = 0 is u(0,t) = 0. The exact 
solution to this problem is especially simple. A shock front 
propagates with constant speed ac, beginning at the closed end, 
and the fluid behind the front is at rest with uniform pressure 
p-^_ = apcv Q . Meanings of the symbols introduced are: 

p = fluid density, 

c = acoustic velocity. 

Factor a is given by 

a = [1 + (ce o /2v Q ) 2 ]' 2 - (ce o /2v Q ). 

In the region ahead of the shock front (x > act) the variables 
p, u, and e maintain their initial values. 

2.5 Governing Equations . Details concerning the governing 
equations are given in Appendix A. Equations are stated in 
forms suitable for any number (i.e., one, two, or three) of 
spatial dimensions. Four separate formulations are developed 
with the principal dependent variable being particle displacement, 
fluid pressure, velocity potential, and displacement potential 
for the respective cases. The capability of each to deal with 
the one - dimens ional water-hammer problem is discussed separately 
below. 

2.4 Displacement Formulation . In this case it is convenient to 
replace the displacement vector by r = p5 and, instead of the 
dilatation e, a density - weighted condensation s = -pe is intro- 
duced. The working equations (Eqs. All, A12, and A10) consist 
of an equation expressing r in terms of the gradient of p, one 
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giving s as the negative of the divergence of r, and the pair of 
algebraic relations for calculating p from s (the bilinear con- 
stitutive law). It is readily seen that the initial conditions 
of the water-hammer problem allow determination of the initial 
values of r and r. The first may be deduced through spatial 
integration of the constant initial dilatation e Q and the second 
is known from the given initial velocity -v . Given these 
required initial values and the boundary condition at x = 0, 
the working equations suffice to solve the problem. 

The displacement vector is continuous, but both velocity 
and pressure are discontinuous at the shock front. The pressure 
discontinuity is not representable by the shape functions of the 
finite element method. Accordingly there must be a finite 
length interval over which both the pressure rise and particle 
deceleration take place. The necessity for this compromise must 
be considered a defect (but not a disqualifying one) of the dis- 
placement formulation. A further disadvantage, not shared by 
any of the other formulations is that the principal dependent 
variable is a vector, not a scalar. In the axisymmetric appli- 
cations planned this will double the order and bandwidth of the 
stiffness matrix:, resulting in a great increase in requirements 
for computer storage and processing time. 

2,5 Pressure Formulation . In Appendix A the pressure formulation 
is stated in terms of a second order equation (A 13) expressing 
s as the Laplacian of the dynamic component of pressure, together 
with the bilinear constitutive relation (A 10) . The solution 
technique involves stepwise time integration for s, alternating 



7 



with use of the constitutive relation to find new values of p. 

In view of the role played by s , it is slightly misleading to 
call this the pressure formulation. Indeed, in the absence of 
cavitation it is advantageous to use the constitutive relation 
to eliminate s in favor of p alone. In cavitated regions, how- 
ever, retaining s makes the strategem of a nonzero bulk modulus 
unnecessary. Despite this advantage, both versions fail. 

The reason this technique fails when applied to the water- 
hammer problem is readily apparent. The variable s (and also p) 
has a step discontinuity at the shock front. (Correspondingly, 
the second derivative s has a dipole singularity.) We know 
that the height of the step depends on the initial (negative) 
value of s and on the fluid velocity. The velocity, however, 
is neither explicitly nor implicitly represented in the equations. 
There seems to be no further reason to consider the pressure 
formulation for cavitation problems. 

The bilinear fluid model with a nonvanishing bulk modulus 
in the tension region was initially introduced with the expec- 
tation that the pressure formulation would work. Since the 
sole reason for including this complication has vanished, we 
shall henceforth exclude negative pressures and correspondingly 
choose 3 = 0 in (A2) , (A4) , and (A10) . 

2.6 Velocity Potential Formulation . Using 4 to represent the 
velocity potential, the governing equations express <j> as the 
negative of the dynamic pressure (A14) , s as the negative of the 
Laplacian of 6 (A15) , and p in terms of s (A10) . In this 
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instance the velocity potential is continuous, but its first 



derivatives are not. It is possible to obtain a moderately 
satisfactory solution from the discretized equations. 

2.7 Displacement Potential Formulation . Using ip to represent 
the displacement potential, the governing equations express ip 

as the negative of the dynamic pressure (A16) , s as the negative 
of the Laplacian of ip (A17) , and p in terms of s (A10) . All of 
the needed variables appear, explicitly or implicitly, in this 
formulation. Moreover, ip and its first derivatives are contin- 
uous. The step discontinuities in p, s, and u are manifested 
as discontinuities in the second derivatives of ip. 

2.8 Formulation Selection . Among the formulations considered, 
the pressure -based one is discarded as unworkable and the velo- 
city potential is rejected as inferior to the displacement poten- 
tial on the basis of continuity. In the remainder of this report, 
discussion of application details will be confined to the ip 
formulation because it is novel. The displacement formulation 
was developed in parallel and also tested on the Bleich -Sandler 
example (Ref. 4) . 

2.9 Discretized Equations and Solution Process. The process of 
forming discretized finite element equations from the correspond- 
ing partial differential equation is well known (e.g., see Ref. 7) 
and will not be repeated here. We note that, based on prior 
experience with wave propagation problems, linear shape functions 
were chosen. It was found that a lumped "mass" matrix was easier 
to use and gave better performance than its consistent counterpart. 
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Details concerning the formulation of initial conditions 
and the radiation boundary condition are given in Appendix B. 
The considerations used for selecting a time integration algor- 
ithm and the means for introducing needed artificial damping 
are discussed in Appendix C. 



3. The Bleich-Sandler Example 

3.1 Statement of the Problem . In Ref. 4 Bleich and Sandler, 
using a bilinear fluid model, study cavitation phenomena during 
one-dimensional wave propagation. They use the method of 
characteristics and introduce additional relations to connect 
state variables on opposite sides of a shock front. 

The numerical example they give concerns "the response of a 
horizontal layer of mass on the surface of a half-space of fluid, 
Fig. 1. A plane pressure wave with a sudden rise and an expon- 
ential decay moves toward the surface, reaching the mass at the 



Atmospheric 





Pressure History 



Figure 1. Particulars 

time t = 0. The system 
pressure, all particles 
shock. The analysis is 



of Bleich-Sandler example (from Ref. 4). 

is subject to gravity and atmospheric 
being at rest prior to arrival of the 
based on the degenerate model with Q = 0." 
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For calculation and presentation of results Bleich and Sandler 
use the time constant of pressure wave decay (- 1 ms) as the 
time unit. The acoustic velocity c is given the convenient 
value unity by choosing unit length to be the distance traveled 
by the pressure wave in unit time (- 56 in.). 

3.2 Comparisons with Bleich-Sandler Results . Bleich and 
Sandler present two figures summarizing their solution. First 
of these traces the time history of the cavitated region in the 
x - t plane. Their figure is reproduced in Figs. 2 and 3 below 
with superposed points obtained by present analyses. For the 
finite element analyses the discretized region extends from 
x = 0 to a radiation boundary at x = 4. Results shown in Fig. 2 




Figure 2. Bleich-Sandler example: time-history of the 
cavitated region. Discrete points found by 
finite element method using the displacement 
formulation . 
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were obtained from the displacement formulation using a node 
spacing (element length) Ax = .04 and a time step h = .01. In 
the absence of damping it was found that moving "islands" of 
positive pressure appeared within the cavitated region. This 
behavior is henceforth called "frothing." It was found that 
damping with r) = .16 was sufficient to supress frothing and 
produce a smooth variation of the condensation s within the 
cavitated region. 

In Fig. 3 the superposed points were obtained using the 
displacement potential formulation. For these results: Ax = .01 




Figure 3. B le ich- Sandler example: time-history of the cavitated 
region. Discrete points found by finite element 
method using the displacement potential formulation. 
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and h = .0025. To suppress frothing a value n = .0025 was 
found to be sufficient. 

A second Bleich-Sandler figure shows the time history of 
the velocity of the surface mass layer. This is reproduced as 
Figs. '4 and 5. Points obtained from present analysis using the 
displacement formulation are superposed on Fig. 4 and points 




Figure 4. Bleich-Sandler example: nondimens ional upward velocity 
of surface mass. Discrete points found by finite 
element method using the displacement formulation. 



from the displacement potential formulation on Fig. 5, The 
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Figure 5. Bleich-Sandler example: nondimensional upward 

velocity of surface mass. Discrete points found 
by finite element method using the displacement 
potential formulation. 



parameters used were those given for Figs. 2 and 3, 
It is believed that the agreement demonstrated 
substantiates the adequacy of displacement and disp 
potential formulations for one - dimens ional analyses 



respectively . 
in Figs . 2-5 
lacement 
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4. Effect of Cavitation on Resilient Attenuator Performance 

We consider here the effect of cavitation on the perform- 
ance of two selected attenuators. One is a GFC with an initial 
thickness of 10 in. under hydrostatic pressure. The second is 

an REL with modulus = 1000 psi. and an initial thickness 

* 

L q = 10 in. when loaded only by atmospheric pressure . The 
shock loading consists of a pressure step of amount p g followed 
by exponential decay with time constant 25 ms. The hydrostatic 
pressure is p^. Results are given as the quotient of the maxi- 
mum dynamic pressure increment Ap max by the peak shock pres- 
sure p . 
r s 

If cavitation effects are ignored the response may be found 
by integrating numerically a first order ordinary differential 
equation. Details are omitted here. 

The response considering cavitation has been determined 
using finite element modelling based on the displacement poten- 
tial formulation. Results are summarized in Table 1. 

Table 1. Effect of Cavitation on Attenuator Response 



I 


Values of Ap /p„ 

r max *s 


Ph 


Ps 


j 

GFC i REL 


psia 


psia 


* 

N.C . 


w . c . * * ! N.C. 1 

I 


w.c. 


30 


1000 


0.21 


0.69 | 0.71 


0.76 


280 


750 


0.44 

i 


0.45 0.82 


0.82 



*N.C. = no cavitation **W.C. = with cavitation 



It is postulated that the relation between gage pressure 
p and thickness L is p = C, (L /L-l) . 

o § ^ 
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Results in Table 1 show that cavitation has little effect 
on the performance of the REL considered at either hydrostatic 
pressure. This is also true for the GFC at p^ = 280 psi., but 
there is severe performance degradation at p^ = 50 psia. Note, 
however, that performance remains better than that of the REL. 

If the hydrostatic pressure is increased significantly above 
280 psia., maintaining the relation + P s = const., cavitation 
will not occur. 



5. Conclusions 

Both the displacement formulation and the displacement 
potential formulation have been shown to produce acceptable 
results when applied to the Bleich -Sandler example. It is 
anticipated that either formulation will provide a workable 
basis for solving three-dimensional axisymmetric problems. 

The scalar displacement potential will lead to a much smaller 
computer storage requirement and processing time, but it may 
not be easy to fit it to the framework of an existing program 
such as NASTRAN, MARC, or NONSAP. No insurmountable difficulty 
is anticipated in using the displacement formulation with one 
of these programs. 



16 



Appendix A. Governing Equations 



Equations are derived here in a form independent of the 
number of spatial dimensions. 

Newton's Second Law: p<$_ = -Vp + f. (Al) 

In the above: 

p = fluid density; 

<5 = particle displacement vector; 

V = gradient operator; 
p = fluid absolute pressure; 
f = body force per unit volume. 

Note that the underline is used to denote a vector quantity. 
Differentiation with respect to time is denoted by a superior 
dot and the convective contributions to the material derivative 
are neglected. 

2 

Bilinear Constitutive Law: p = -c pe, e <_ 0; 

C A 2 ) 

p = -Bcpe, e>0. 



Here : 

c = acoustic velocity in fluid; 

e = dilatation. 

2 

Note that c p is the bulk modulus of the fluid. For the bilinear 
fluid model 6 is chosen as positive and small compared with unity. 
The limiting condition of zero pressure in the cavitated region 
corresponds to g = 0 . 

Geometric identity: e = V’<5_ , (A3) 

where the dot denotes the scalar product. 
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It is possible to choose a single dependent variable 
such as p and, through suitably chosen manipulations, demon- 
strate that p obeys the wave equation in the uncavitated fluid 
and a modified form with 3c in place of c in the cavitated 
region(s) . Thus 

2 „ 2 ^ n 

p = cVp, p >_ 0 ; 

.. 7 7 7 (A4) 

p = 3 c V p, p < 0. 



A more enlightening approach which focuses attention on 
the sequential steps in time integration of the governing 
equations uses auxiliary dependent variables and a set of 
three equations. For this purpose we first define some addi- 
tional dependent variables and then summarize four separate 
formulations . 

Definitions . In our applications the body force f appearing 
in (Al) may be expressed as 

£ = Vp h , (A5) 

where p^ is the hydrostatic component of fluid pressure. 

It is useful to introduce two density weighted variables: 

r = p6_ , (A6) 

s = -pe. (A7) 

We also introduce two similarly weighted potential 
functions : 



= r, (AS) 

Vip = r. (A9) 

Henceforth we omit explicit reference to the density factor and 
refer to r as displacement, s as condensation (Lamb's usage. 
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